R Introduction for Biology

Dealing with data

The tidyverse is an eco-system of packages that provides a consistent, intuitive system for data manipulation and visualisation in R.

A spreadsheet that we would normally analyse in Excel can can be imported into R by using a function from the readr package.

• read_csv
• read_tsv

The readxl package can also import files in Excel’s own format

• read_xls
• read_xlsx

We will use the file tcga_clinical_CLEANED.tsv (11MB) as an example. This file contains information about biological samples that were included as part of the TCGA (The Cancer Genome Atlas) project - all tumour types.

If you do not have this file, it can be downloaded with the following code. N.B. the code will first check if the file is present or not.

if(!file.exists("tcga_clinical_CLEANED.tsv")) download.file(url = "https://sbc.shef.ac.uk/r4biol/tcga_clinical_CLEANED.tsv", 
              destfile = "tcga_clinical_CLEANED.tsv")

As the file name ends in .tsv the function we want is probably read_tsv

library(readr)
clinical <- read_tsv("tcga_clinical_CLEANED.tsv")

If you got an error could not find function "read_tsv" at the previous step, it will because you don’t have the readr package installed. In which you will need to run the command:-

install.packages("readr")
library(readr)

This will install the readr package. You will then need to run the read_tsv line of code from above.

If you get really stuck reading data into R, you can use the Import Dataset option from the File menu which will allow you to choose the parameters to read the data interactively

The Environment panel of RStudio (top-right) should show that an object called clinical has been created. This means that we can start doing some analysis on this dataset. The choice to call the object clinical was ours. We could have used any name instead of clinical but we chose something vaguely informative and memorable.

The dimensions of the object should should 7706 obs. of 420 variables. This means the object we have created contains 7706 rows and 420 columns. Each row is a different observation; in this case a different biological sample. Each column records the value of a different variable. In R’s terminology, we have just created a tibble which is a special case of something called a data.frame. As we will see the the object can contain either numbers or text in each column.

RStudio has a browsing functionality available that we can activate using the View function

View(clinical)

Observations that are missing are indicated by an NA. It may look at first glance that the table mostly contains missing data. This is because some columns are only recorded for particular tumour types and not all tumour types recorded the same information.

Data Exploration

We will now start to explore the data with a couple of different questions in mind:-

• What are the differences in male/female split between different tumour types?
• Is there a difference in diagnosis age in different tumour types?

The tasks required to do this will be greatly helped by the dplyr package.

library(dplyr)

As with the readr package above, you will need to have installed the dplyr package beforehand.

install.packages("dplyr")

Choosing what columns to analyse

The select function allows us to narrow down the number of variables we are interested in from 420. The first argument is always the name of the data frame. There are numerous different ways of specifying which column(s) you want, including typing the names of the columns of interest in exactly the same way that they appear in the data. Let’s assume we already know the names of columns containing tumour type and gender.

The dataset utilizes a binary classification of gender as ‘male’ or ‘female’. It is important to note that this categorization may not fully encompass the diverse range of gender identities recognized today.

##Note that the spelling of "tumor" has to exactly match that found in the data
select(clinical, 
       tumor_tissue_site,
       gender
       )

We can reasonably guess that the Age information would be contained in a column that has “age” somewhere in the name. However, the column is not just called “age” or “Age”

select(clinical, age)
select(clinical, Age)

Without manually going through the columns, there are a few “helper” functions that we can employ

select(clinical, contains("age"))

select(clinical, contains("age_"))

select(clinical, starts_with("age"))

Up to now we have not changed the underlying dataset. select is showing what the dataset looks like after applying the specified subset. If we want to make permanent changes we can create a variable

analysis_data <- select(clinical,
                        bcr_patient_barcode,
                        tumor_tissue_site,
                        gender,
                        age_at_diagnosis)

Restricting rows

The select function only performs the very specific task of letting you choose what columns you want to analyse. After using select, our dataset analysis_data still has all 2496 rows.

The function to choose or restrict to the rows we might be interested in is called filter. We have to write a short R command to choose the rows.

e.g. if we want only the male samples we use the following code. Notice that two “=” signs are required. If you try and use the function with a single “=” R will print a helpful hint.

filter(analysis_data, gender == "MALE")

and the females

filter(analysis_data, gender == "FEMALE")

An equivalent, but somewhat unnecessary for this example, statement would be to ask for rows where the gender is not equal to MALE

filter(analysis_data, gender != "MALE")

Identifying male or females based on the gender column is simplified by there only being two possible entries in this column; MALE or FEMALE. In a real life situation, especially if data have been entered by hand, there could be inconsistencies that require action. See an example at the end of the tutorial.

We can also restrict the data based on numeric columns by using >, < etc.

filter(analysis_data, age_at_diagnosis > 80)

But what if we want males with Brain tumours? dplyr allows us to combine more than one condition if we separate them with a ,. In computing this is known as an “and” statement and only rows where both statements are two will be shown. The line could be extended to include more than two tests if we wanted to.

filter(analysis_data, gender == "MALE",tumor_tissue_site == "Brain")

How about brain or lung tumours? Using a | symbol instead of a , allows for either of two (or more) conditions to be TRUE.

filter(analysis_data, tumor_tissue_site == "Brain" | tumor_tissue_site == "Lung")

To answer the question of how many males / females have a certain tumour type we could now use R statements such as:-

filter(analysis_data, gender == "MALE",tumor_tissue_site == "Brain")

filter(analysis_data, gender == "FEMALE",tumor_tissue_site == "Brain")

and make a note of the number of observations included in the resulting data frame. However, there is much more flexible way of summarising data in this manner.

Summarising

Although useful for data exploration, it would clearly be inefficient to get gender/tumour type counts in this way as we would have to repeat for all combinations of tumour type and gender. The function count can now give us exactly what we want. The output is given as a tibble, so we could use some of the functions that we have learnt about so far (select, filter…) to further manipulate. e.g. obtain the counts for just Brain/

count(analysis_data, 
      tumor_tissue_site,
      gender)

The count function is useful for tabulating the number of observations, but for other summary statistics a more general sumamrise function can be used. This can be used in conjunction with basic summary functions supported by base R. A summary statistic being something that can be applied to a series of numbers and produce a single number as a result. e.g. the average, minimum, maximum etc.

summarise(analysis_data, 
          Average = mean(age_at_diagnosis),
          min = min(age_at_diagnosis),
          max = max(age_at_diagnosis))

However, we have a problem due to missing values. If R sees and missing values in a column it will report the mean, minimum or maximum of that column as a missing value. Although this default behaviour can be changed, before proceeding we could also choose to remove any missing observations from the data. These are represented by a NA value, which is a special value and not a character label.

filter(analysis_data, is.na(age_at_diagnosis) | is.na(tumor_tissue_site))

We want the opposite of the above; where the age of diagnosis and tumour site is not missing.

analysis_data <- filter(analysis_data, !is.na(age_at_diagnosis), !is.na(tumor_tissue_site))

The summary will now work as expected.

summarise(analysis_data, 
          Average = mean(age_at_diagnosis),
          min = min(age_at_diagnosis),
          max = max(age_at_diagnosis))

This might not be what we want in all circumstances, as the statistics can also be calculated on a per-tumour site basis using dplyrs group_by function.

data_grouped <- group_by(analysis_data, tumor_tissue_site)
summarise(data_grouped,Average = mean(age_at_diagnosis),
          min = min(age_at_diagnosis),
          max = max(age_at_diagnosis) )

Sorting (arranging)

We have previously used filter to restrict the rows that we are interested in. Rather than just analysing the male or female patients (for example), we might also want the rows in our table to be ordered according to the gender column.

arrange(analysis_data, gender)

We can also arrange by columns containing numeric values in either ascending (the default) or descending order.

arrange(analysis_data, age_at_diagnosis)

## Use a descending order
arrange(analysis_data, desc(age_at_diagnosis))

Like how sorting works in Excel, we can also use mutliple columns for sorting. e.g. if we want ordering by diagnosis age for each tumour type separately.

arrange(analysis_data, tumor_tissue_site, age_at_diagnosis)

Workflows and “piping”

So far we have used several operations in isolation. However, the real joy (?) of dplyr is how different operations can be chained together. Lets say we just wanted female tumours.

filter(analysis_data, gender == "FEMALE")

Our next step could be to remove the gender column since it is somewhat redundant.

analysis_data2 <- filter(analysis_data, gender == "FEMALE")
select(analysis_data2, tumor_tissue_site, age_at_diagnosis)

## or 
## select(analysis_data2, -gender)

The code would quickly get cumbersome if we wanted to include additional steps such as removing NA values. An alternative approach called “piping” is recommended and activated by adding %>% at the end of a line. This tells R to use the output of the current line as the first argument on the next line. In this current example it means we don’t need to specify which data frame that select uses as input - it will use the data frame created by the filter in the previous line. The code written using %>% is more concise.

filter(analysis_data, gender == "FEMALE") %>% ## and then...
  select(tumor_tissue_site, age_at_diagnosis) ## %>% and then...

The %>% operation becomes available when you load dplyr. If you wish to use piping outside of dplyr there is also a “base” equivalent |> that doesn’t require any libraries to be loaded

filter(analysis_data, gender == "FEMALE") |> ## and then...
  select(tumor_tissue_site, age_at_diagnosis) ## |> and then...

We recently created a summary table for each tumour type giving the average, minimum and maximum of diagnosis age. This can be replicated using %>% and an extra sorting step added to the end.

group_by(analysis_data, tumor_tissue_site) %>% 
summarise(Average = mean(age_at_diagnosis),
          min = min(age_at_diagnosis),
          max = max(age_at_diagnosis)) %>% 
  arrange(Average)

Overview of plotting

Our recommending way of creating plots in RStudio is to use the ggplot2 package - especially as it interacts well with dplyr and other tidyverse packages.

library(ggplot2)

A couple of useful references are given here-

• ggplot2 reference guide
• Flowchart for deciding on what graph type to use

The general principle of creating a plot is the same regardless of what kind of plot we want to make

• specify the data frame containing the data we want to plot
• specify which columns in that data frame we want to use for various aesthetic aspects of the plot
• define the type of plot we want
• apply any additional format changes

A bar plot would be a natural choice for showing the counts of male / female samples. The geom_bar plot will automatically count how many occurrences there are for each value.

ggplot(analysis_data, aes(x = gender)) + geom_bar()

Numerical data can be visualised using a density plot or histogram. The density is automatically calculated and displayed on the y-axis.

ggplot(analysis_data, aes(x = age_at_diagnosis)) + geom_density()

In order to compare the age distributions of different tumour types we can also imagine this being displayed as a series of boxplots with

• the age variable on the y-axis
• the type of tumour on the x-axis

this can be translated into ggplot2 language as follows -

ggplot(analysis_data, aes(x = tumor_tissue_site, y = age_at_diagnosis)) + geom_boxplot()

A disadvantage of the boxplot is that it only gives a very crude summary of the data. It can be misleading when applied to data with few observations and is often preferable to add individual data points

ggplot(analysis_data, aes(x = tumor_tissue_site, y = age_at_diagnosis)) + geom_boxplot() + geom_jitter(width=0.1)

Adding some colour to the plot can be achieved by adding a fill aesthetic and specifying what column to map the colours too. A colour palette is automatically chosen, but can be changed afterwards if we wish.

ggplot(analysis_data, aes(x = tumor_tissue_site, y = age_at_diagnosis, fill = tumor_tissue_site)) + geom_boxplot() + geom_jitter(width=0.1)

Adding the fill aesthetic for the density plot can be used to show a separate curve for each tumour type.

## alpha of 0.5 used to make the curves transparent
ggplot(analysis_data, aes(x = age_at_diagnosis, fill = tumor_tissue_site)) + geom_density(alpha=0.5)

Another useful technique for splitting the plots based on a variable is to use the facet_wrap function that will give a grid of plots. For instance we can show male/female counts for each tumour type separately.

ggplot(analysis_data, aes(x = gender,fill=gender)) + geom_bar() + facet_wrap(~tumor_tissue_site)

By combining all the techniques we have seen we can compare the diagnosis age between males and females; separately for each tumour type.

ggplot(analysis_data, aes(x =gender, y = age_at_diagnosis, fill = gender)) + geom_boxplot() + geom_jitter(width=0.1) + facet_wrap(~tumor_tissue_site)

Challenges of “messy” data

Real-life data are often less straightforward to deal with than the “cleaned” dataset presented here. Despite the many high-throughput technologies that are used for scientific investigation, there is inevitably a spreadsheet(s) needed to describe the experimental setup and this is typically entered manually.

So-called “Data Wrangling” is a crucial and time-consuming part of the analysis process taking 80% of analysis time by some estimates. Hadley Wickham, Chief Scientist at Posit and lead author of ggplot2 likens tidy and messy data to Leo Tolstoy’s quote about families:-

Happy families are all alike; every unhappy family is unhappy in its own way

Like families, tidy datasets are all alike but every messy dataset is messy in its own way.

A comprehensive guide to the issues surrounding data entry via spreadsheets, and how to avoid them, is given by Data Carpentry.

• Data Carpentry Spreadsheets lesson

However, for public data that we have no control over we often have no choice but to clean the data ourselves. We have intentionally created an alternative dataset with a few intentional issues to illustrate the cleaning process. The following code will download the data if you do not have it already.

if(!file.exists("tcga_clinical_MESSY.tsv")) download.file(url = "https://sbc.shef.ac.uk/r4biol/tcga_clinical_MESSY.tsv", 
              destfile = "tcga_clinical_MESSY.tsv")

messy <- read_tsv("tcga_clinical_MESSY.tsv")
messy

Whitespace

“whitespace” is the addition of a blank character or space to the beginning or end of text. Traditionally it is a problem because it will create extra categories in your data. e.g. MALE and MALE. The messy dataset that you have just imported includes some whitespace in the tumor_tissue_site column. However, the read_tsv function automatically ignores whitespace values as the trim_ws argument of read_tsv is set to TRUE (see the help page ?read_tsv).

messy_ws <- read_tsv("tcga_clinical_MESSY.tsv", 
                     trim_ws = FALSE)
messy_ws

count(messy_ws,tumor_tissue_site)

The resulting data frame now contains two apparently identical categories for Bladder. However, with the use of the nchar function, which counts the number of characters, we can see that extra spaces must be included.

count(messy_ws,tumor_tissue_site) %>% 
  mutate("Length_of_Label" = nchar(tumor_tissue_site))

Clearly we could have used the default settings for read_tsv and the problem would not have occurred. Otherwise, it is useful to know about the stringr package that contains many useful functions for cleaning character data.

• The stringr package

For the example of removing whitespace we can use the str_trim function combined with a mutate. This will replace all the whitespace in the tumor_tissue_site column and overwrite the column. If we repeat a count afterwards we see only the unique entries that we expect.

library(stringr)
mutate(messy_ws, tumor_tissue_site = str_trim(tumor_tissue_site)) %>%
  count(tumor_tissue_site)

Inconsistent coding of variables

Unfortunately the tumor_tissue_site column is not the only one with issue that need fixing with these data. If, as before, we try and plot the number of males/females in the dataset we get a surprise.

ggplot(messy, aes(x = gender)) + geom_bar()

There is no differentiation between female and FEMALE or male and MALE. Whilst we can intuitively decide that these represent the same value, they do not get automatically combined in R. The consequence being that attempts to identify all the male patients will require some careful coding. The example used previously will now no longer identify all the correct patients.

filter(messy, gender == "MALE")

One solution when filtering would be to add different criteria to account for the different capitalisation

filter(messy, gender == "MALE" | gender == "male")

However, since we know that the error is due to inconsistent use of upper/lowercase we can use the str_to_upper function in stringr to convert all values to uppercase. Or indeed we could convert all to lowercase using str_to_lower if we prefered.

messy %>% 
  mutate(gender = str_to_upper(gender)) %>% 
  ggplot(aes(x = gender)) + geom_bar()

We would probably want to also make the change permanent by creating a new variable

cleaned <- read_tsv("tcga_clinical_MESSY.tsv") %>% 
    mutate(gender = str_to_upper(gender)) 

A more generic approach would be to use the forcats package to replace all occurrences of male with MALE and the same for females.

• The forcats package

The package allows us to “recode” entries in a column that contains a factor. i.e. categorical.

library(forcats)
mutate(messy, gender = forcats::fct_recode(gender,"MALE"="male"),
       gender = forcats::fct_recode(gender,"FEMALE"="female")) %>% 
    ggplot(aes(x = gender)) + geom_bar()

The approach is more flexible than merely changing the case, as it could also replace other values such as “m” or “f” if they existed.

## Just example code if we wanted to replace "m" with "MALE"
messy <- messy %>%
  mutate(gender = fct_recode(gender,
                            "MALE" = "male", 
                            "MALE" = "m" 
  ))

Different means of representing missing values

There are many different strategies for representing missing data. The read_tsv function should automatically detect any NA values in the source dataset and treat them appropriately. However, in our messy dataset we also have NULL values (as seen by making a count of the values in age_at_diagnosis).

count(messy, age_at_diagnosis) %>% arrange(desc(n))

Because the NULL value is present in the age_at_diagnosis column, R will treat. the entire column as containing characters. Therefore we cannot use the kind of plots we would expect with numeric data

ggplot(messy, aes(x = age_at_diagnosis)) + geom_histogram()

Likewise we can’t calculate numeric summaries; although R will attempt to and create a data frame of NA values rather than giving an error.

  group_by(messy, tumor_tissue_site) %>% 
  summarise(Mean_Diagnosis_Age = mean(age_at_diagnosis,na.rm=TRUE))

Warning: There were 19 warnings in `summarise()`.
The first warning was:
ℹ In argument: `Mean_Diagnosis_Age = mean(age_at_diagnosis, na.rm = TRUE)`.
ℹ In group 1: `tumor_tissue_site = "Bladder"`.
Caused by warning in `mean.default()`:
! argument is not numeric or logical: returning NA
ℹ Run ]8;;ide:run:dplyr::last_dplyr_warnings()dplyr::last_dplyr_warnings()]8;; to see the 18 remaining warnings.

The read_tsv function has the ability to replace NA values when the data are imported. Specifically, the na argument can be used to define what values are being used to represent missing.

read_tsv("tcga_clinical_MESSY.tsv", na = c("NULL","NA")) %>% 
  ggplot(aes(x = age_at_diagnosis)) + geom_histogram()

read_tsv("tcga_clinical_MESSY.tsv", na = c("NULL","NA")) %>% 
  group_by(tumor_tissue_site) %>% 
  summarise(Mean_Diagnosis_Age = mean(age_at_diagnosis,na.rm=TRUE))

Including units in the column

The final column in this example contains height information (where available) for our patients. Clearly it is important to know what units this is recorded in, but placing the units inside the entries creates issues as we can’t treat the data as numbers.

arrange(messy,height_at_diagnosis)

The stringr package can be used again, and this time a function called str_remove_all which removes all occurrences of a particular string. In particular we want to remove cm. We will need an additional step to convert the column into numeric values

messy %>% 
  mutate(height_at_diagnosis=str_remove_all(height_at_diagnosis, "cm")) %>% 
  mutate(height_at_diagnosis = as.numeric(height_at_diagnosis)) %>% 
  arrange(height_at_diagnosis)

There is usually more than one way of completing a task in R. In this instance, we could also use the str_sub function in stringr to extract a “substring” from each entry in the column. The argument end=-3 specifies that the extraction should end three characters before the end of each string.

messy %>% 
  mutate(height_at_diagnosis=str_sub(height_at_diagnosis, end=-3)) %>% 
    mutate(height_at_diagnosis = as.numeric(height_at_diagnosis)) %>% 
  arrange(height_at_diagnosis)

Final code to clean the data

For reference, here is the final code chunk that can be used to clean the data.

cleaned <- read_tsv("tcga_clinical_MESSY.tsv", na = c("NULL","NA")) %>% 
  mutate(messy, gender = forcats::fct_recode(gender,"MALE"="male"),
       gender = forcats::fct_recode(gender,"FEMALE"="female")) %>% 
    mutate(height_at_diagnosis=str_sub(height_at_diagnosis, end=-3)) %>% 
    mutate(height_at_diagnosis = as.numeric(height_at_diagnosis))